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Abstract 

Integrating and differentiating matrices allow the numerical integration 
and differentiation of functions whose values are known at points of a 
discrete grid. Previous derivations of these matrices have been restricted to 
one-dimensional grids or to rectangular grids with uniform spacing in at least 
one direction. The present work develops integrating and differentiating 
matrices for grids with non-uniform spacing in both directions. The use of 
these matrices as operators to reformulate boundary value problems on 
rectangular domains as matrix problems for a finite-dimensional solution 
vector is considered. The method requires non-uniform grids which include 
"near-boundary" points. An eigenvalue problem for the transverse vibrations 
of a simply supported rectangular plate is solved to illustrate the method. 
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1 • Introduction 


Rotating beam configurations have traditionally been used to study the 
vibrations and aeroelastic stability of rotating structures such as helicopter 
rotor blades and propeller blades. More recently, models involving elastic 
plates have been proposed to include the effects of spanwise variations in 
material properties. The fourth-order boundary value problems associated with 
both the beam and plate models do not, in general, have useful closed form 
solutions. Consequently, most theoretical work on these problems has been 
asymptotic or numerical in nature. 

In one approach to the numerical solution of these problems, harmonic time 
dependence is assumed to reduce the governing partial differential equation to 
a differential equation in space variables which includes an eigenvalue. For 
beam models, this is an ordinary differential equation. The fundamental 
derivative which represents beam curvature may now be taken as a new dependent 
variable, and the eigenvalue problem for the beam can be reformulated as an 
integro-differential equation (White & Malatino, 1975; Kvaternik, White, & 
Kaza, 1978; Lakin, 1982). This equation may be conveniently expressed using 
integral, differential, and boundary evaluation operators. The operator 
equation for the continuous solution may further be converted to a matrix 
operator equation for a finite-dimensional solution vector by evaluating the 
continuous equation at a finite set of discrete grid points which span the 
interval of interest. A key question is now the manner in which the matrix 
operators are approximated. 

For beam models, one method for approximating the integral and 
differential operators involves the use of integrating matrices (Vakhitov, 
1966; Hunter, 1970; Lakin, 1979) and differentiating matrices (Hunter & 
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Jainchell, 1969; Lakin, 1985). In the simplest terms, these matrices provide, 
respectively, a means of numerically integrating and differentiating a 
function whose values are known at a finite set of discrete grid points. A 
key property of both integrating and differentiating matrices is that their 
derivation requires only knowledge of the grid points, and no information is 
needed about the function to be numerically integrated or differentiated. In 
the case of a beam model with its single space variable, this property allows 
the integrating and differentiating matrices based on one-dimensional grids to 
be used directly as approximations for the integral and differential operators 
in the matrix operator form of the eigenvalue problem. The result of this 
approximation is a striaghtf orward matrix eigenvalue problem which can be 
solved by standard methods. This approach has proved capable of efficiently 
handling a wide variety of beam problems including beams with concentrated 
masses, follower forces, and point loadings (Lakin, 1982). 

For vibration and buckling problems which involve two-dimensional elastic 
plates, removal of the time dependence from the original boundary value 
problem yields an eigenvalue problem which continues to be governed by a 
partial differential equation. By analogy with the one-dimensional case, it 
would seem desirable to reformulate this eigenvalue problem as a matrix 
integro-partial differential equation for a finite-dimensional solution vector 
on a two-dimensional grid of discrete points. Integrating and differentiating 
matrices based on two-dimensional grids could then be used to approximate the 
respective operators resulting, again, in a standard matrix eigenvalue 
problem. 

The present work will explore the potential of this approach by 
considering an eigenvalue problem associated with the transverse vibration of 
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a simply supported rectangular plate. This problem consists of the biharraonic 
eigenvalue equation 


a 2 J2. ~ o 

( — j + — o) u(x,y) = X u(x,y) (1.1) 

8x 9y 

for 0 x A , ^ _< y _< B, and t ^ ie boundary conditions 

u(0,y) = u (0,y) = 0, 
u(A,y) = u xx (A,y) = 0, 

( 1 . 2 ) 

u(x,0) = u yy (x,0) = 0, 
u(x,B) = u (x,B) = 0. 

yy 

In section 2, equation (1.1) will be reformulated as an integro-partial 
differential equation consistent with the form of the boundary conditions 
(1.2). Because (1.2) involves conditions on u(x,y) itself at all four 
boundaries, the present approach retains u itself as the dependent 
variable. Conversion to a matrix eigenvalue problem will now require the 
derivation of appropriate integrating and differentiating matrices based on a 
two-dimensional rectangular grid of discrete points. 

One type of integrating matrix for a function of two variables has been 
previously derived by Lakin (1981). This matrix may be used in two- 

dimensional problems whose reformulation is possible using an integrating 
matrix alone, e.g., the plate analogue of a beam with cantilevered boundary 
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conditions. Unfortunately, this matrix is not suitable for the present 
purposes as its derivation requires that the spacing of the grid points be 
uniform in at least one direction. The boundary conditions (1.2) lead to a 
reformulation which will require the use of differentiating matrices to 
approximate partial derivatives with respect to both x and y. To preserve 
accuracy, the computational grid must now include points "close" to all four 
boundaries (Lakin, 1985) giving non-uniform grid spacing in both directions. 

Generalized integrating matrices for functions of two variables whose 
values are known on non-uniform rectangular grids are derived in the 
Appendix. Differentiating matrices which approximate partial derivatives on 
non-uniform rectangular grids are also derived, as are matrices which evaluate 
quantities at boundary grid points. The only restriction in this derivation 
is that the grid sub-units [ Xj < x < x j+J , y k < y < y^] should be 
rectangles. 

In section 3, appropriate integrating and differentiating matrices are 
used to approximate the reformulated eigenvalue problem by a matrix eigenvalue 
problem involving a "stacked" column vector. Numerical calculations are 
presented for a rectangle with A = 2 and B = 1 on a 7-by-7 grid 
(including near boundary points). Despite the coarseness of this grid, good 
agreement with the exact eigenvalues of (1.1) and (1.2) is obtained. 


2. Reformulation of the Eigenvalue Problem 

Before considering the two-dimensional problem (1.1) and (1.2), it is 
useful to briefly consider the reformulation process for the one-dimensional 
problem of a cantilevered beam. For the beam, integrating the fourth-order 



ordinary differential equation twice with respect to the space variable x 
(say) gives an integral equation for the curvature w"(x) and allows the 
boundary conditions at the free end to be explicitly invoked. Boundary 
conditions at the free end now enter explicitly through relation of the 
fundamental derivative w" to w' and the original dependent variable w. 
Thus, the appropriate reformulation makes use of all boundary conditions. 

The initial steps required to obtain an appropriate reformulation of the 
two-dimensional problem (1.1) and (1.2) are a straightforward generalization 
of the one-dimensional procedure. Equation (1.1) may be integrated twice with 
respect to x from x to A making use of the boundary condition u xx (A,y) 
= 0. This result may then be integrated two additional times with respect 
to y from y to B making use of three other boundary conditions at x = A 
and y = B , i . e . , 

u(A,y) = u(x,B) = u yy (x,B) = 0. (2.1) 

These steps give the equation 
B 

/ (5 “ y>[ u xx < x »S) + (A “ x)u xxx (A,£)]d5 

y 

A 

+ / (n - x)[u yy (n,y) + (B - x)u yyy (n,B)]dn 

(2.2) 

+ 2(A - x) [u (A,y) - u (A,B)] + 2(B - y)[u (x,B) - u (A,B)] 

xx y y 

2 B A 

+ 2 (A - x) (B - y)u (A,B) + 2u(x,y) = X / (g - y) J (n - x)u(n,C)dn d£ . 

x y 



The corner consistency condition u(A,B) = 0 has also been explicitly used in 
obtaining (2.2). 

Equation (2.2) contains three different types of second partial 
derivatives. A question to be answered is therefore whether, as in the single 
variable beam case, the boundary conditions at x = 0 and y = 0 can be 
accounted for by designating one type of second derivative as a fundamental 
derivative and new dependent variable. The form of the boundary conditions is 
critical in deciding this question. 

The second partial derivatives u xx and Uyy are related to the 
variable u through the relations 

x 

/ (n ~ x)u xx (n,y)dn = u(x,y) - u(0,y) - xu x (0,y) (2.3) 

and 

y 

/ (C - y)u yy (x,5)d£ = u(x,y) - u(x,0) - yu y (x,0). (2.4) 

Thus, for either u xx or Uyy alone to be a viable candidate for the 
fundamental derivative role, the boundary conditions (1.2) would have to 
specify that the respective first partial derivative vanish for x or y 
equal zero. 

In contrast to (2.3) and (2.4) the boundary conditions in (1.2) at 
x — 0 and y — 0 insure that integrated terms vanish when the mixed second 
derivative u is related to u through the expression 

y x 

0 0 U xy^ n, ^ dn = u (x,y) - u(x,0) - u(0,y) + u(0,0). 


(2.5) 
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This is reraeaiscent of the single variable case. Unfortunately, the higher 
derivative conditions at x = 0 and y « 0 are not on u xy itself but on 
u xx and Uy y • The mixed derivative can be directly related to u xx and 
u yy through an appropriate quadrature and partial differentiation. However, 
the relationship is such that neither of the two remaining higher derivative 
boundary conditions can be satisfied directly. Consequently, the mixed 
derivative alone is also not suitable as a new dependent variable. 

Two possibilities now remain. First, the boundary value problem (1.1) and 

(1.2) can be reformulated in terras of three fundamental derivatives, i.e., all 
three types of second partial derivatives can be designated as separate 
dependent variables. The original variable u could then be related to the 
new variables through (2.5) satisfying two of the four conditions at x = 0 
and y = 0. However, while attractive in principle, this option will be 
impractical in practice as it leads to matrix eigenvalue problems involving 
large matrices. For example, in the case of the relatively coarse 7 -by- 7 
rectangular grid used in the present calculations for the sample problem, 147- 
by-147 matrices would be required if the three derivatives are retained as 
distinct dependent variables. 

A second possibility is to retain the original variable u(x,y) in the 
reformulation of (1.1) and (1.2). This choice requires that two additional 
partial derivative operations be approximated. It also requires that equation 

(2.2) be modified to explicitly take account of the four boundary conditions 
at x = 0 and y — 0. However, when the resulting reformulation is 
approximated by a matrix eigenvalue problem for a f inite-dimensional solution 
vector, considerably smaller matrices are needed than if fundamental 
derivatives were used. For example, in calculations on the 7-by-7 
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rectangalar grid, 49-by-49 matrices are sufficient for the basic development, 
and the final eigenvalue problem involves a 25-by-25 matrix. 

To obtain the desired reformulation, equation (2.2) is first evaluated 
at x = 0, at y = 0, and at the corner point x = 0 and y = 0. The three 
relations which result now explicitly use all conditions at the boundaries 
(0,y), (x,0), and satisfy the corner consistency conditions at (0,0). In 

particular: 

A B 

/ I i[u yy (n,y) + (B - y)u yyy (n,B)]dq + / A(? - y) u xxx (A,5)d5 


+ 2A[u x (A,y) - u x (A,B)] - 2(B - y)[u y (0,B) - u y (A,B)] (2.6) 

2 B A 

+ 2A(B - y)u (A,B) = X / (5 “ y) / Tiu(n,5)dq d£ , 
y y 0 


B A 

/ S[ u xx (x, 5 ) - (A - x)u xxx (A,5)]d? + J B(q - x)u (n,B)dq 
0 v j y y 


+ 2 (A - x)[u x (A,0) - u x (A,B)] + 2B[u y (x,B) - u (A,B)] (2.7) 

2 B A 

+ 2B(A - x)u (A,B) » X J £ f (i) - x)u(q,0 dp d£ , 

7 Ox 


B A 

/ A Cu (A, 5 )d§ + J Bqu (q,B)dn 
0 0 yyy 


+ 2A[u x (A,0) - u x (A,B)] + 2B[u y (0,B) - u (A,B)] (2.8) 

B A 

+ 2ABu (A,B) = X / £ / r\u(n,Odn d£ . 

y 0 0 
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Subtracting (2.6) and (2.7) from (2.2) and adding (2.8) now gives 

!/ 5 + y / |[*u CA,E) - u xxx ( x ,5)]d S 

(0 y ) 

+ lo n + *1 }^ yu yyy <n,B) “ u yyy (n,y) ^ dn 

+ 2x[u x (A,0) - u x (A,y)] + 2y[u y (0,B) - u y (x,B)] 

(2.9) 

+ 2xyu (A,B) + 2u(x,y) 

2 ( y A B x 

“ x ) X J / 5u(n,5)dri dJ; + yj / nu(r),5)dn d£ 

( 0 x y 0 

BA y x \ 

+ xy / / u(n,Odn d S + / / n5u(n,?)dn d£ ( . 

y x 0 0 ) 

This equation is the reformulation required for appoximation of the eigenvalue 
problem (1.1) and (1.2) using integrating and differentiating matrices on the 
rectangular grid. 

3 . Approximation By A Matrix Eigenvalue Problem 

The first step in approximating equation (2.9) by a matrix eigenvalue 
problem for a finite-dimensional solution vector is to discretize this 
equation on a rectangular grid G of discrete points. This will allow the 
integrating, differentiating, and boundary matrices on this two-dimensional 
grid to be used as operators to approximate the corresponding operators in the 
continuous equation. 
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The two-dimensional grid G may be formed from the cross-product of 
appropriate one-dimensional grids in the x- and y-directions . In 

particular, let G x be the one-dimensional grid of N discrete points 

0 = Xj < x 2 ••• < x N = A (3.1) 

which discretizes the interval 0 < x < A, and let G y be the one-dimensional 
grid of M points 

0 = y l < y 2 < < y M = B (3.2) 

which discretizes the interval 0 < y < B. Neither G„ nor G„ need to have 
uniform spacing. Indeed, in actual implementation, to preserve accuracy it 
will be necessary to choose the spacings - x^ , x^ - x^_^ , , and 

y M “ relatively small as the formulation will involve differentiating 

matrices (Lakin, 1985). The two-dimensional grid G for the continuous 
region of the boundary value problem (1.1) and (1.2) may now be taken as the 
set of NM discrete points 


G = {(^jyj) :x i e G x »yj€Gy} • (3.3) 

Thus, the subunits of the grid G are rectangles which need not have equal 
areas . 

The straightforward format for displaying values of the solution u(x,y) 
at the discrete grid points of G is an N-by-M matrix [U] with elements 
U ij = uUjL.yj). Unfortunately, this format is unsuitable for the desired 
reduction to a matrix eigenvalue problem. Rather, it is convenient to arrange 
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the NM values of u on the grid as a "stacked" NM-by-1 column vector. In 
particular, the finite-dimensional solution vector {u} is taken to have 
elements 

u, = u(x ,y.) with k = N(j - 1) + i (3.4) 

K i J 

where i « 1,...,N and j = 1,...,M. It should be noted that this format 
for {u} is x-oriented, i.e., the M groups of N consecutive elements in 
the stacked vector give values of u for a fixed value of y in Gy while 
x varies in G x . Because of this orientation, as indicated in Appendix A, 
matrices which approximate integrals and derivatives with respect to y will 
require extra operations in their construction. A general flow chart for the 
construction of both x- and y-operation matrices is given in Figure 1. 

Once equation (2.9) has been discretized, integrating, differentiating, 
and boundary evaluation matrices on the grid G may be used to obtain a 

matrix eigenvalue problem which provides the required approximations. In 

particular, the eigenvalue problem for (2.9) may be written as 

[G] {u } = X 2 [H]{u} (3.5) 

where [G] = [G x ] + [G 2 ] + 2 [G 3 ], 

[GJ = ( [ JRY ] [Y] + [Y ] [ JTY ] ) ( [X] [BA] [DX3] - [DX2 ] ) , (3.6) 

[G 2 ] = ( [ JRX] [X] + [X] [ JTX] ) ( [Y] [ BB ] [DY3] - [DY2]), (3.7) 


and 
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[G 3 1 = [ X ] ( [ BAO ] - [ BA ] ) [ DX ] + [Y]([BBO] - [BB])[DY] 


+ [X] [Y] [BA] [BB] [DX] [DY] + [I]. 


(3.8) 


In (3.6) through (3.8), [JRX] and [JTX] approximate x— integrals from 0 
to x and from x to A, repectively, on the two-dimensional grid G 
while [DX] , [DX2], and [DX3] approximate first, second, and third partial 
derivatives with respect to x on G. Similarly, [ JRY ] and [JTY] 
approximate y-integrals from 0 to y and from y to B, respectively, 
on G while [DY], [DY2], and [DY3] approximate first, second, and third 
partial derivatives with respect to y. The matrices [BA] and [BB] 
evaluate quantities at boundary points with x = A and y = B, respectively, 
while [BAO] and [BBO] give values at the boundary points (A,0) and 

(0,B). The derivation of these integrating, partial differentiating, and 
boundary matrices on the grid G is discussed in Appendix A. [I] is the 
NM-by-NM identity matrix. The NM-by-NM matrices fX] and [Y] are 

^^■^S^nal matrices such that if position k in the stacked solution vector 
corresponds to the point (x^y^) of G, then X kk = x< and Y kk = y y 
The matrix [H] in (3.4) may be written as the sum of four matrices 
[H] = [Hj] + [H 2 ] + [H^] + [H^] with 


[Hj] 

= [X] [JRY] [JTX] [Y], 

(3.9) 

[h 2 ] 

= [Y] [JTY] [JRX] [X], 

(3.10) 


[H 3 ] = [X] [Y] [JTY ] [JTX] , 


(3.11) 
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and 


[H 4 ] = [JRY ] [ JRX] [X] [Y] . 


(3.12) 


To enhance accuracy in calculation of the lower eigenvalues, it is convenient 
to further re-write (3.5) as 


[A] |u} = to {u} (3.13) 

with 

[A] = [G] -1 [H] and w = l/\ 2 . (3.14) 

To test the accuracy of this matrix eigenvalue approximation to the 
continuous problem (1*1) and (1*2), equation (3*13) was solved on a rectangle 
with A = 2 and B = 1. G x was taken to be the seven point grid consisting 
of the two end points x^ = 0 and x-j = 2, two near-boundary points X2 - 
0.0001 and x & =1.9999, and three interior points X3 = 0.5, = 1.0, and 

X5 = 1.5. Gy was taken to be the seven point grid consisting of the two 
endpoints y^ = 0 and y 7 = 1.0, two near-boundary points y2 = 0.00005 
and y^ = 0.99995, and three interior points y3 = 0.25, y^ = 0.5, and 
y 5 = 0.75. The implementation of (3.6) through (3.12) on this grid thus 
involves 49-by-49 matrices. However, the size of the matrix [A] in the 
matrix eigenvalue problem (3.13) can be reduced by noting that u = 0 at all 
boundary points. Rows and columns of [A] corresponding to boundary points 
may thus be deleted leading to a 25— by— 25 matrix. 

Exact solution of the boundary value problem (1.1) and (1.2) are of the 


form 


u (x,y) = c sin sin 
pq pq A B 


(3.15) 
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where p and q are integers. The corresponding eigenvalues are 


X 


pq 



(3.16) 


For the present test case with A - 2 and B = 1, the exact values of the two 
smallest eigenvalues are 


X u = 12.337 and X 21 = 19.739. (3.17) 

Approximate values for these eigenvalues were obtained by solving the matrix 

eigenvalue problem (3.13) on the grid G. Differentiating matrices were based 

on fourth degree polynomials while integrating marices were based on fifth 
degree polynomials. This even/odd degree scheme allows grid points to be 
centered as much as possible within the sliding subgrids on which the 

differentiating and integrating matrices are based. The computations give the 
values 

A n = 12.553 and X 21 = 17.635. (3.18) 

It must be remarked that G in this test problem Is a relatively coarse grid 
which, when boundary and near-boundary points are ommitted, has only a total 
of nine points in the interior of the two-dimensional region. Values of these 
approximations could be improved through the inclusion of additional interior 
points. 
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4* Concluding Remarks 

The present work has examined an extension of integrating and 
differentiating matrix methodology to partial differential equations involving 
two space variables. Matrices which approximate integrals and derivatives on 
one-dimensional grids are used as a starting point to develop matrices which 
approximate integrals and partial derivatives on two-dimensional rectangular 
grids. The method requires that the original boundary value problem be 
reformulated to take account of all boundary conditions. Integrating, 
differentiating, and boundary matrices may then be used as operators to 
approximate the boundary value problem by a standard matrix problem for a 
stacked, finite-dimensional solution vector. The inclusion of near boundary 
points in the grid helps to prevent the degradation of accuracy at boundaries 
associated with differentiating matrices. 

While only two-dimensional rectangular domains have been explicitly 
considered in the present work, a further generalization of the method to 
three-dimensional domains is relatively straightforward. This is due to the 
use of a stacked column vector format for the solution vector which allows 
matrices on the higher dimensional grid to be obtained from matrices on the 
underlying one-dimensional grids. The primary requirement in going to three 
space dimensions is the use of appropriate change matrices analogous to the 
matrices [CXY] and [CYX] of Appendix A. These matrices will shuffle the 
order of the stack to orient it with respect to a given variable and then 
restore the original orientation after an operation with respect to that 
variable has been approximated. 
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Appendix A: Partial Differentiating, Integrating, and Boundary Matrices 

For Rectangular Grids 

Let G be the rectangular grid in (3.3) formed from the two one- 
dimensional grids G x and Gy in (3.1) and (3.2). Let {u} be the 
NM-by-1 "stacked" column vector defined in (3.4) which gives values of 
u(x,y) at the points of G. Then, for the present boundary value problem on 
a rectangle, it is necessary to derive two types of integrating matrices and 
three partial differentiating matrices for each of the two space variables 
x and y . 

Consider first the matrices which approximate operations with respect to 

x. The "0-to-x" integrating matrix [JRX] on G is an NM-by-NM matrix 

such that the NM-by-1 column vector [JRX]{u} contains approximate values 

x 

of the integral / u(r\,y)dn at the points of G. In particular, 

0 

X 2 A X 2 

[JRX] {u} ~ (0,/ u(n ,y. )dp ,••*,/ u(n,y- )dn,0,/ u(n ,y 9 )dn ,• • • , 

0 0 0 

(Al) 

A x 2 A T 

/ u(n,y 9 )dn,... ,0,/ u(n,y )dn,*‘*,J u(n,y )dn) . 

0 0 n 0 


Similarly, the "x-to-A" integrating matrix [JTX] leads to approximations 

A 

at points of G of the integral f u(n,y)dn so that 

x 


AAA 
[jtx] { u} ~ (/ u(n ,y. )dn >•••>/ u(n,y. )dn,0,/ u(n,y 7 )dn, 

" 0 Vi 0 

A A 

/ u(n,y M )dn,***,J 

0 


,o,. 


(A2) 


u(n,y M )dn,0) T . 



-17- 


The x-partial derivative matrix [DX] is an NM-by-NM matrix such that 
[DX] {u} contains approximate values of 3u/3x at points of G. Thus, 


[ DX ]{ U } _ (|£ (x N , 7l >,•••, <Vy M >. •••»!£ (x N’ y M^ T * 


The matrices [DX2] and [DX3] which lead to approximations of second and 
third partial derivatives of u(x,y) with respect to x at points of G 
have similar definitions. 

The stacked column vector {u} has been constructed so as to be x- 
oriented. It thus consists of M segments which contain N elements 
apiece. In each of these segments, x varies through G x for a constant 
value of y in G y . The matrices in (Al) to (A3) approximate x-operations 
for fixed values of y. Consequently, use of the present stacked column 
vector format will allow x— operation matrices on the rectangular grid G to 
be constructed from the corresponding N-by-N matrices on the one-dimensional 
grid G x . This construction is most easily accomplished through definition of 
a "diagonalizing" mapping from the set of N-by-N matrices to the set of NM- 
by-NM matrices. 

Let p and q be integers, and let [A] be a p-by-p matrix. The 
diagonalizing mapping Diag(p,q, [A] ) then assigns to [A] the pq-by-pq 
matrix [B] obtained by placing q matrices [A] along the diagonal of 
[B] and taking all other elements of [B] to be zero, i.e., 

[A], 0 

[B] = Diag(p,q , [A] ) = [A] . # . (A4) 

_° *[A]_ 
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The x-operation matrices on G may be formed by applying this mapping with 
p = N and q = M to appropriate matrices on the grid G. 


Consider first the construction of 

[JRX]. 

Let 

f(x) 

be a 

function whose 

values are known at the points of G, 

and let 

{f} 

be 

the 

N-by-1 

column 

vector which contains these values. 

Further, 

let 

[jrx] 

be an 

N-by-N 

integrating matrix which approximates 

integrals 

of 

f(x) 

from 

0 to 

x on 


G x so that 


[jrx] {f } 



f(n)dn 


(A5) 


Comparing equation (Al), segraent-by-segment , with (A5) now shows that 


[JRX] = Diag(N,M, [jrx] ). (A6) 

The matrix [jrx] (and hence [JRX]) is not unique, but depends on both the 
number of points included in the sliding subgrids of G x and the manner in 
which f(x) is approximated on these subgrids. The present work uses the 
integrating matrix [jrx] for one-dimensional non-uniform grids developed by 
Lakin (1979). 

To obtain the second required integrating matrix [JTX] on the 

rectangular grid G, let [jtx] be an integrating matrix on G x such that 

[jtx] {f } ~ 

Then, a segment-by-segraent comparison of (A2) and (A7) shows that 


A 

/ f(n)dn 




(A7) 



[JTX] = Diag(N ,M, [ j tx] ) . 


(A8) 


Dif f erentiating matrices on one-dimensional non-uniform grids have been 
derived by Lakin (1985). Let the matrices which approximate first, second, 
and third derivatives of f(x) at points of G x be denoted, respectively, 
by [dx], [dx2], and [dx3]. As is the case with integrating matrices, these 
differentiating matrices are not unique. Further, in the usual case where the 
sliding subgrids contain fewer than all N points of G x , [dx2] and [dx3] 
cannot be obtained by simply squaring or cubing the matrix [dx]. Rather, 
these matrices must be obtained directly from approximations to the second and 
third derivatives of f(x) on the sliding subgrids. 

The three matrices [DX], [DX2], and [DX3] which approximate partial 
derivatives with respect to x of u(x,y) at points of G may now be 
constructed from [dx], [dx2], and [dx3] using the diagonalizing mapping. 

In particular, 

[DX] = Diag(N,M, [dx] ) , 

[DX2] = Diag(N,M, [dx2] ) , 

and (A9) 

[DX3] = Diag(N,M, [dx3 ] ) • 

For consistency, the differentiating matrices used in the present work were 
based on the same subgrid approximation scheme as was used for the integrating 
matrices on G x . 

Consider next derivation of the matrices which approximate operations with 
respect to y on the rectangular grid G. The required integrating matrices 
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for this variable are the "O-to-y" integrating matrix [JRY] which 

y 

approximates the integral J u(x,5)d? at points of G, and the "y-to-B" 

0 

B 

integrating matrix [JTY] which approximates the integral J u(x,£)d£. The 

y 

product of these NM-by-NM matrices with {u} are the NM-by-1 column 

vectors 


y 2 y 2 

[JRY]{u} ~ (0, •••,(),/ u(x ,5)d5,**«,J u(x 

0 0 

B B 

/ u(x 1 ,5)d?,»»»,/ u(x ,5)d^) T 

0 0 


(A10) 


and 


B B 

t JTY] {u} ~ (f u(x ,C)d £,••*,/ u(x ,£)dg,*”, 

0 0 

B B 

/ u(x. ,5)d? ,•••,/ u(x ,S)d£,0,...,0) T . 

y M-l y M-l 


(All) 


Let g(y) be a function whose values are known at the points of the one- 
dimensional grid Gy, and let {g} be the M-by-1 column vector which 
contains these values. Further, let [jry] and [jty] be integrating 
matrices on Gy such that 


[jry] {g} 


( y i 

J 

, 0 


g(5)d£ 


and 


[jty] {g} ~ 



g(?)d5 


(A12) 


Because {u} , as defined, is x-oriented, the integrating matrices [JRY] 
and [JTY] on G cannot be formed using the single mapping (A4) on the 
corresponding matrices for the one-dimensional grid G y . As [JRY] and 
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[JTY] approximate integrals with respect to y for fixed values of x, two 
additional mappings which convert from x- to y-orientation and from y- 
to x-orientation will also be required. 

If the values of u(x,y) at the NM points of the rectangular grid G 
are written as an NM-by-1 stacked column vector {v} whose format is y- 
oriented, then the elements of {v} are 

v^ = u(x^,yj) where k = M(i - 1) + j. (A13) 

Thus, the vector {v} consists of N segments which contain M elements 
apiece. In each segment, y varies in Gy for a fixed value of x in G x . 
If {v} and not {u} had been chosen as the format for the solution vector, 
then [JRY] and [JTY] could be formed from [jry] and [jty] directly 

using the diagonalizing mapping (A4) with p = M and q = N. 

An x-oriented vector {u} may be associated with its corresponding 
y-oriented vector {v} through a mapping Cxy from the set of NM-by-1 
column vectors into itself so that Cxy({u}) = {v}. In symbolic terms, if the 
values of u(x,y) on the rectangular grid G are arranged in an N-by-M 
array, then the effect of applying the mapping Cxy is to produce an 
M-by-N array which is the transpose of the original. For the present 

purposes, the mapping Cxy may be carried out by multiplying {u} by an NM- 
by-NM matrix [CXY] so that 


[CXY ] {u} = {v} . 


(A14) 


The matrix [CXY] may be written as a stack of N, M-by-NM matrices 
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(A15) 


If le^J is the j-th unit vector in N-dimensional real space, i.e., a row 
vector with a one in the j-th position and zeros in the other N-l 
positions, then each matrix in (A15) can be written in the form 

0 

le.J 

J • . . (A16) 

V 

The matrices in (A15) thus have the row vector [e^ J along the diagonal and 
zeros elsewhere. 

The mapping Cxy is one“to — one and hence invertable. Let the inverse 

mapping be denoted by Cyx so that Cyx({v}) = {u}. Thus, when applied to 
a y-oriented vector, Cyx restores the standard x-orientation. This mapping 
may be carried out by multiplying {v} by the NM-by-NM matrix [CYX] so 
that 

[CYX] {v } = {u}. ( A17 ) 

Because of the inverse relationship of Cxy and Cyx, [CYX] = [CXY]" 1 . 

Having defined the mappings which change the format of a stacked column 
vector from x- to y-orientation and back again, the procedure which uses 
the M-by-M integrating matrix [ j xry ] on Gy to produce the the vector 

[ JRY] {u} defined in (A10) may now be given. The vector {u} is first 
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multiplied by [CXY] to produce a y-oriented format. It is then multiplied 
by the matrix Diag(M,N, [ j ry] ) to give a vector which consists of M 
segments, each of which approximates the integral of u(x,y) from 0 to y 
for a fixed value of x. Finally, the original x-orientation is restored 
through multiplication by [CYX] . This process implies 

[JRY] = [CYX] Diag[M,N, [ jry] ) [CXY] . (A18) 

The matrix [JTY] may likewise be formed in this manner from the integrating 
matrix [jty] on Gy. Hence, 

[JTY] = [CYX] Diag(M,N, [jty]) [CXY ] . (A19) 


The NM-by-NM matrix [DY], which approximates partial derivatives with 
respect to y on the rectangular grid G, is defined by 


[DY H u l = Of- (x N»yi>»' 


8u 


»3y (X 1 ’ 




(x n» y m (A20) 


Let [dy] be an M-by-M differentiating matrix on the one-dimensional grid 
Gy such that 

[dy] {g} ~ {g'} • (A21) 


Then, [DY] can be formed from [dy] through the relation 


[DY] = [CYX] Diag (M , N , [ dy ] ) [ CXY ] . 


(A22) 
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Similarly, let [dy2] and [dy3] be matrices which approximate second and 
third derivatives of g(y) on G y . Then, replacing [dy] in (A22) by 
[dy2] or [dy3] leads to the matrices [DY2] and [DY3] , respectively, 
which approximate second and third partial derivatives of u(x,y) with 
respect to y on G. 

Matrices which evaluate quantities at the boundaries x = A and y = B 
and at the corner points (A,0) and (0,B) are the final items needed to 
construct the matrices [G] and [H] in the matrix reformulation of (1.1) 
and (1.2). The NM-by-NM boundary matrix [BA] is such that [BA]{u} gives 
values of u(A,y) at points of the rectangular grid G. In particular, 


[BA] {u} = (u(A,y 1 ),»*«,u(A,y 1 ),u(A,y 2 ),»»» , 


u (A,y 2 ) > , **>u(A,B),»»» , u(A,B)) T . 


(A23) 


[BA] may be written as a stack of M, N-by-NM matrices. If [baj] is the 
j-th matrix in this stack (j = 1,...,M), then the element in the i-th 

row (i — 1 , . . . ,N) and k— th column of [baj] is unity if k = jN and zero 
otherwise. Similarly, the NM-by-NM matrix [BB] is such that [BB][u] 

gives values of u(x,B) at points of G, i.e.. 


[BB] { u} = (u(x 1 ,B),*.»,u(x N ,B),*..,u(x 1 ,B),»»*,u(x N ,B)) T . (A24) 


[BB] may also be written as a stack of M, N-by-NM matrices. 


[BB] = 


[bb] 

[bb] 


(A25) 
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However, each of the M matrices in the stack (A25) is identical. The right- 
hand block of N columns of [bb] is simply the N— by-N identity matrix. 
All other elements of [bb] are zero. 

The matrices [BAO] and [BBO] evaluate quantities at the corner 
points (A,0) and (0,B), respectively. Both of these NM-by-NM matrices 
consist of a single non-zero column which contains all ones. For [BAO], 
the N-th column is non-zero. The (NM-N+l)-st column is non-zero in the 
case of [BBO], 
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x-operation matrices y-operation matrices 

[JRX] [JTX] [JRY] [JTY] 

[DX] [DX2] [DX3] [DY] CDY2] [DY3] 


Figure 1. 


Flowchart illustrating the formation of matrices which 
approximate integrals and partial derivatives with respect to 
x and y on the rectangular grid G. 
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